C................. AURORA.FOR ...............................
C... This is the main program for the 2-stream model for
C... calculating the auroral electron spectrum. Revised September
C... 1987 by Phil Richards. It reads the control parameters from
C... the command file, sets up the altitude grid, sets up the
C... neutral atmosphere by calling MSIS or reading from file and
C... finally calls the subroutine PE2S to calculate the auroral
C... spectrum.

      SUBROUTINE AURORA(TIME,      !universal time in hours, in
     >	                 ZWR,      !altitude for writing spectrum
     >	              AURFAC,      !auroral factor in/out
     >               DTMAX_S)      !maximum time step,out

      PARAMETER (IZDIM=401)
      INTEGER AURFLAG
      REAL MASS(3),ZWR
      LOGICAL DEBUG
      DOUBLE PRECISION DALT,GLN,GL,GLATM
      DIMENSION XN(3,IZDIM),SIGEX(3),SIGION(3),XN0(3),Z(IZDIM)
     > ,ZTE(IZDIM),XNE(IZDIM),BG(IZDIM),DZ(IZDIM),D(19),T(2),TN(IZDIM)
      DATA MASS,IXS/16,32,28,1/, ISWFLX/-1/
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,UTSEC,GL(401)
      DATA ITCONS,  TE ,XNE0,EVMIN , AVMU  ,STEP,JUBND,DIPANG,ZWRITE
     >   /    0 , 3000, 6E5  , 0.0, 0.577 , -9 , -9  , 1.57 ,300  /
      DATA AVMU/0.577/,DEBUG/.FALSE./

      !...get the auroral factors
      CALL GET_AURFAC(TIME,AURFAC,AURFLAG,ECHAR_S,SFAC_S,TENFLX_S,
     >  DTMAX_S) 

      !.. Find maximum electron energy (EVMAX) for  auroral calculation
      CALL MAX_ENERGY(AURFAC,SFAC_S,ECHAR_S,TENFLX_S,EVMAX)

      !...check the AURFAC and AURFLAG to determine to continue or return.
      !...if AURFAC is too small or AURFLAG is 0,return.Otherwise,continue.
      IF(AURFLAG.EQ.0.OR.AURFAC.LT.0.01) RETURN

      !...set up the altitute grids
      CALL SETUP_ALTGRID(IZDIM,BG,Z,DZ,JUBND,STEP)
      !... Call the MSIS model for the neutral atmosphere
      GLATM=GL(1)
      DO 30 J=1,JUBND
        DALT=Z(J) !.. double precision altitude for AMBS
        CALL AMBS(1,DALT,GLATM,D,T,CHIF,SAT,GLON,GLATD)
        XN(1,J)=D(2)  !.. [O]
        XN(2,J)=D(4)  !.. [O2]
        XN(3,J)=D(3)  !.. [N2]
        TN(J)=T(2)    !.. Tn
        !...... thermal electron densities
        XNE(J)=XNE0
        IF(Z(J).LE.110.0) XNE(J)=XNE0/10.0**((110.-Z(J))/15.)
        IF(Z(J).GT.300.0) XNE(J)=XNE(J-1)*EXP(-DZ(J)/300.0)
        !... The Te profile approaches TN at low altitudes and TE at high
        ZTE(J)=(TN(J)*100.0+(Z(J)-100.0)*TE)/Z(J)
        IF(ZTE(J).LT.TN(J)) ZTE(J)=TN(J)
 30   CONTINUE
      
      !...... write out the neutral atmosphere
      IF(DEBUG) THEN 
        WRITE(6,391)
 391    FORMAT(/2X,' Altitude varaiation of densities and'
     >   ,X,'temperatures'/)
        DO 43 J=1,JUBND
        IF(J.EQ.1.OR.MOD(J,15).EQ.0) WRITE(6,902)
 902    FORMAT(2X,'    J     Z(J)    DZ(J)    Tn     [O]   ',
     >   '  [O2]     [N2]      Te      [e]')
        WRITE(6,90) J,Z(J),DZ(J),TN(J),(XN(IS,J),IS=1,3),ZTE(J),XNE(J)
 43     CONTINUE
      ENDIF

      !..... Call routine to calculate the auroral fluxes
      CALL AUR2STR(IZDIM,Z,XN,JUBND-1,STEP,XNE,ZTE,BG,ZWR,ECHAR_S
     > ,ISWFLX,TN,SFAC_S,EVMAX,EVMIN,AVMU,DIPANG,ITCONS,TENFLX_S,CHI)
      RETURN
 90   FORMAT(2X,I6,3F8.1,1P22E9.1)
 91   FORMAT(I4,F7.2,1P7E10.2,0P22F8.1)
 92   FORMAT(1P22E10.2)
      END


C:::::::::::::::::::::::::GET_AURFAC:::::::::::::::::::::::::::
C......This subroutine is called by AURORA to get the AURFAC.
      SUBROUTINE GET_AURFAC(TIME,   !universal time in hours, in
     >	                  AURFAC, !auroral factor, out
     >                      AURFLAG,!a flag to indicate change of parameters
     >                      ECHAR_S,!ECHAR scalar
     >                      SFAC_S, !SFAC  scalar
     >                      TENFLX_S, !TENFLX scalar
     >                      DTMAX_S)  !DTMAX scalar

      
      INTEGER JCALL               !counts number of calls
      INTEGER INDEX               !index of time in last call
      DATA JCALL/0/               !initial value=0
      DATA INDEX/1/               !initial value=1
      INTEGER INDEX_NEW           !index of time in current call
      DATA INDEX_NEW/1/           !initial value=1

      !data structures for calling READ_AUR_DATA()
      INTEGER ALTDIM              !... Array dimensions 
      PARAMETER (ALTDIM=999)
      REAL UTORLT                 !... UT OR LT
      REAL UT(ALTDIM)             !... Universial time
      REAL TENFLX(ALTDIM)         !... Engery flux
      REAL ECHAR(ALTDIM)          !... Characteristic electron energy
      REAL SFAC(ALTDIM)           !... Shape factor of Gaussian or Maxwellian.
      REAL DTMAX(ALTDIM)          !... Maximum time step in seconds
      !data to be return to AURORA()
      REAL ECHAR_S                !... ECHAR scalar
      REAL SFAC_S                 !... SFAC  scalar
      REAL TENFLX_S               !... TENFLX scalar
      REAL DTMAX_S                !... DTMAX scalar
      INTEGER AURFLAG
      AURFLAG=0 
                  
      JCALL=JCALL+1  !.. number of calls to this routine

      !..If this is the first time, read the auroral data file
      IF(JCALL.EQ.1) THEN
         AURFAC=0.0
         CALL READ_AUR_DATA(ALTDIM,MAXDIM,UTORLT,UT,TENFLX,ECHAR,
     >      SFAC,DTMAX)

      !..also do error checking 
         CALL ERR_CHECK(MAXDIM,ZLBDY,
     >        UTORLT,UT,TENFLX,ECHAR,SFAC,DTMAX)
      ENDIF

      !..find index of input time in time array
      CALL FIND_INDEX(TIME,UT,INDEX,MAXDIM,INDEX_NEW)

      !.. If TIME not found in arrays set AURFAC to zero and return
      IF (INDEX_NEW.EQ.-1000) THEN
          AURFAC=0.0
          RETURN
      ENDIF
       
      !... copy auroral parameters to scalars to return their values
      ECHAR_S=ECHAR(INDEX_NEW)
      SFAC_S=SFAC(INDEX_NEW)
      TENFLX_S=TENFLX(INDEX_NEW)*6.25E11        ! convert to eV
      DTMAX_S=DTMAX(INDEX_NEW)

      !.. Check if index same as previous one. If yes,return previous AURFAC      
      IF (INDEX_NEW.EQ.INDEX.AND.JCALL.GT.1) RETURN  
           
      !.. Current index is not same as previous. Check to see if there 
      !.. is an aurora	    
      IF(TENFLX(INDEX_NEW).LT.0.01.OR.ECHAR(INDEX_NEW).LT.1.0) THEN
         AURFAC=0.0
      ELSE
        AURFAC=1.0
        !..Check parameters to see if any have changed enough.
        !..If yes, return AURFAC. If no, reevaluate auroral
        IF(ABS(TENFLX(INDEX_NEW)-TENFLX(INDEX)).GT.0.01) AURFLAG=1
        IF(ABS(ECHAR(INDEX_NEW)-ECHAR(INDEX)).GT.0.01) AURFLAG=1
        IF(ABS(SFAC(INDEX_NEW)-SFAC(INDEX)).GT.0.0001) AURFLAG=1
      ENDIF        

      INDEX=INDEX_NEW       !..update index for next time
      
      RETURN                                              
      END

C:::::::::::::::::::::::::	SETUP_ALTGRID:::::::::::::::::::::::::::
C......This subroutine is called by AURORA to set up the altitude grid.
      SUBROUTINE SETUP_ALTGRID(IZDIM, !.. array dimensions   
     >                         BG,    !.. grid scaling factor
     >                         Z,     !.. auroral altitudes
     >                         DZ,    !.. altitude steps
     >                         JUBND, !.. index of upper boundary altitude
     >                         STEP)  !.. basic altitude step
      DIMENSION Z(IZDIM),BG(IZDIM),DZ(IZDIM)
      !....... initialize first altitude point to lower boundary
      STEP = 0.25
      ZMAX = 800.00 
      Z(1)= 80.00

      !........ set up the altitude grid .......................
      !.. Cubic grid spacing. BG is needed to calculate derivatives
      !.. DZ= altitude step to next point at altitude Z. By putting
      !.. BG(J)=1.0 the altitude steps are constant. This happens if
      !.. too large a STEP is given

      DO J=1,IZDIM-1
      BG(J)=1E6/(Z(J)*Z(J)*Z(J))
      IF(STEP.GE.0.5) BG(J)=1.0
      DZ(J)=STEP/BG(J)
      Z(J+1)=Z(J)+DZ(J)
      IF(Z(J).LE.ZMAX) JUBND=J
      IF(Z(J).GT.ZMAX+100.) GO TO 933
      
      ENDDO
 933  CONTINUE
      RETURN 
      END


C:::::::::::::::::::::::::FIND_INDEX:::::::::::::::::::::::::::
C......This subroutine is called by AURORA to search the index
C......of the input time in time array.

      SUBROUTINE FIND_INDEX (TIME,     !time to be searched, in
     >                       UT,       !time array,in
     >                       INDEX,    !previous index, in
     >                       MAXDIM,   !max index in time array,in
     >                       INDEX_NEW)!index of input time,out
      

      INTEGER DIMALT
      PARAMETER (DIMALT=999)
      REAL UT(DIMALT) 
      !......do sequential search from last index
      
      DO 300 I=INDEX,MAXDIM
      !......search from previous index to MAXDIM-1
      
         IF(I.LT.MAXDIM) THEN
             IF (TIME.GE.UT(I).AND.TIME.LT.UT(I+1)) THEN
                 INDEX_NEW=I
                 RETURN
             ELSE 
                 GO TO 300
             ENDIF
      !......check last element in time array
         ELSE
             IF(I.EQ.MAXDIM.AND.TIME.EQ.UT(I)) THEN
                 INDEX_NEW=I      
      		   RETURN
             ELSE                 !if input time goes beyond range,
                 INDEX_NEW=-1000  !return a special value to indicate this
                 RETURN
      	   ENDIF	    
         ENDIF
 300  CONTINUE
      END

C:::::::::::::::::::: READ_AUR_DATA :::::::::::::::::::::::::::
C... This subroutine actually reads the AE data from the file
C... Data must either increase monotonically ordecrease monotonically
      SUBROUTINE READ_AUR_DATA(ALTDIM,  !... Array dimension in
     >                        MAXIND,   !... Max index of array out
     >                        UTORLT,   !... 1.0=UT,0.0=LT out
     >                        UT,       !... Universal time out
     >                        TENFLX,   !... Energy flux out
     >                        ECHAR,    !... Characteristic electron energy out
     >                        SFAC,     !... Shape factor out
     >	                    DTMAX)    !... Maximum time step in seconds out
              
      IMPLICIT NONE
      
      INTEGER ALTDIM         !... Array dimensions in GET_AUR_DATA
      INTEGER DIMALT         !... Array dimensions in READ_AUR_DATA
      INTEGER MAXIND         !... Max index of array
      PARAMETER (DIMALT=999)
      INTEGER I              !... loop control variable
      INTEGER FEXISTS        !... Indicates if file exists
      REAL UTORLT            !... UT OR LT
      REAL UT(DIMALT)        !... Universial time
      REAL TENFLX(DIMALT)    !... Engery flux
      REAL ECHAR(DIMALT)     !... Characteristic electron energy
      REAL SFAC(DIMALT)      !... Shape factor of Gaussian or Maxwellian.
      REAL DTMAX(DIMALT)     !....Maximum time step in seconds

      MAXIND=0 0      !... index of last parameter value

      !... Using DIMALT in here because ALTDIM generates a large
      !... executable file for some reason
      IF(DIMALT.NE.ALTDIM) THEN
         WRITE(6,*) ' Incorrect dimension for DIMALT in READ_AUR_DATA'
         CALL RUN_ERROR    !.. print ERROR warning in output files
         STOP
      ENDIF

      !... see if auroral data file exists
      CALL TFILE(14,FEXISTS)
      IF(FEXISTS.EQ.0) RETURN

      !... remove comments before reading data columns
      CALL REM_COMMENTS(14)
   
      !... Read the AE data into  arrays
      DO I=1,DIMALT
        READ(14,*,ERR=96,END=95) UT(I),TENFLX(I),ECHAR(I),
     >    SFAC(I),DTMAX(I)
      ENDDO

 95   CONTINUE

      !.. Make sure all auroral data are read in
      IF(I.GE.DIMALT) THEN
        WRITE(6,199) DIMALT
        CALL RUN_ERROR    !.. print ERROR warning in output files
        STOP
      ENDIF
 199  FORMAT(/'  ** ERROR: Too many times in AURORA parameter file;'
     > ,1X,'must be fewer than ',I5' lines of data in file'/)
 
      MAXIND=I-1   !... Record the maximum index

      UTORLT=1.0 !.. Currently no provision for LT
      !... determine if time is UT or LT
      !IF(ABS(UT(MAXIND)-24.0).LT.0.01) THEN
      !  UT(MAXIND)=24.0
      !  UTORLT=1.0
      !  WRITE(6,*) '  *** Local Time assumed in auroral file'
      !ENDIF

      REWIND(14)
      RETURN

 96   CONTINUE
      WRITE(6,*) '   *** ERROR: in Aurora (unit FOR014) data file'
 100  FORMAT(1X,F8.2,F5.1,F8.0,F8.0,F5.2,F8.2)
      CALL RUN_ERROR    !.. print ERROR warning in output files
      STOP
      END
      
C::::::::::::::::::::::REM_COMMENTS::::::::::::::::::::::::::::::::::::
C.... This routine is used to remove the comment line from data file.
C.... Assume all comment lines starting with '!'
      SUBROUTINE REM_COMMENTS(INUNIT)    !... Input unit in
      CHARACTER*80  ACH                  !... Character array
      
 50   READ(INUNIT,'(A80)',END=60) ACH
      !... if it is comment line, go back to read next line
      DO I=1,80
      if (ACH(I:I).NE.' ') go to 70
      ENDDO
 70   IF (ACH(I:I).EQ.'!') GOTO 50
      !... if it is not comment line,go back one record
      BACKSPACE INUNIT
 60   RETURN
      END


C::::::::::::::::::::ERR_CHECK :::::::::::::::::::::::::::
C... This subroutine actually checks the validity of AE data. 
C... UT must increase monotonically.If it is local time, it
C...    should be less than or equal to 24.00.
C... TENFLX should be within (0.0,100.0).
C... ECHAR  should be within (0,100,000).
C... SFAC   should be within (0.0,10.00).
      SUBROUTINE ERR_CHECK(MAXIND,   !... Max index of array in
     >                     ZLBDY,    !... Lower boundary altitude in
     >                     UTORLT,   !... 1.0=UT,0.0=LT in
     >                     UT,       !... Universal time in
     >                     TENFLX,   !... Energy flux in
     >                     ECHAR,    !... Characteristic electron energy in
     >                     SFAC,     !... Shape factor in
     >                     DTMAX)    !... Maximum time step in seconds	    
      
      INTEGER MAXIND         !... Max index of array
      INTEGER I              !... loop control variable
       
      INTEGER DIMALT         !... Array dimension
      PARAMETER (DIMALT=999)
    
      REAL ZLBDY             !... Lower boundary altitude
      REAL UTORLT            !... Universal time or local time   
      REAL UT(DIMALT)        !... Universal time
      REAL TENFLX(DIMALT)    !... Engery flux
      REAL ECHAR(DIMALT)     !... Characteristic electron energy
      REAL SFAC(DIMALT)      !... Shape factor of Gaussian or Maxwellian.
      REAL DTMAX(DIMALT)     !... Maximum time step in seconds


      !... check UTORLT
      IF (UTORLT.NE.1.0 .AND. UTORLT.NE.0.0) THEN
          WRITE(6,*) 'UTORLT is invalid in aurora data file!'
          CALL RUN_ERROR    !.. print ERROR warning in output files
          STOP
      ENDIF
      !... check UT is monotoniclly increasing or not.
      DO I=1,MAXIND
         IF(I.GT.1.AND.UT(I).LT.UT(I-1)) THEN
            WRITE(6,*)' UT should increase monotonically'
            WRITE(6,*)' in aurora data file'
            CALL RUN_ERROR    !.. print ERROR warning in output files
            STOP
         ENDIF
      !... check UT is less than or equal to 24.00 if it is local time
         IF(UTORLT.EQ.0.0.AND.UT(I).GT.24.00) THEN
            WRITE(6,*)' local time should be less than or equal '
            WRITE(6,*)' to 24.00 in aurora data file'
            CALL RUN_ERROR    !.. print ERROR warning in output files
            STOP
         ENDIF
      !... check TENFLX is within(0.0,100.0) or not
         IF(TENFLX(I).LT.0.0.OR.TENFLX(I).GT.100.0) THEN
            WRITE(6,*) I, 'th TENFLX data wrong in aurora data file!'
            CALL RUN_ERROR    !.. print ERROR warning in output files
            STOP
         ENDIF
      !... check ECHAR is within (0,100000) or not 
         IF(ECHAR(I).LT.0.0.OR.ECHAR(I).GT.100000.0) THEN
            WRITE(6,*) I, 'th ECHAR data wrong in aurora data file!'    
            CALL RUN_ERROR    !.. print ERROR warning in output files
            STOP
         ENDIF
      !... check SFAC is within (0.00,10.00) or not
         IF (SFAC(I).LT.0.00.OR.SFAC(I).GT.10.00) THEN
            WRITE(6,*) I, 'th SFAC data wrong in aurora data file!'
            CALL RUN_ERROR    !.. print ERROR warning in output files
            STOP
         ENDIF
      ENDDO
      RETURN
      END      
C::::::::::::::::::::::::::::: MAX_ENERGY ::::::::::::::::::::::::::::::::::::::::::
C.. Adjust maximum electron energy for  auroral calculation
      SUBROUTINE MAX_ENERGY(AURFAC,SFAC_S,ECHAR_S,TENFLX_S,EVMAX)
      IMPLICIT NONE
      REAL AURFAC,SFAC_S,ECHAR_S,TENFLX_S,EVMAX

      !.. Adjust maximum electron energy for Gaussian auroral calculation
      EVMAX=(3.0*SFAC_S+1)*ECHAR_S
      !.. Adjust maximum electron energy for Maxwellian auroral calculation
      IF(SFAC_S.LT.0.01) EVMAX=NINT(10*ECHAR_S)

      !... make sure energy is not too big
      IF(ECHAR_S.GT.200000) THEN
          WRITE(6,94) ECHAR_S
          CALL RUN_ERROR    !.. print ERROR warning in output files
          STOP
      ENDIF
 94   FORMAT(//' Auroral characteristic energy ',F12.1,
     >    ' eV is too large')

      IF(AURFAC.GT.0.01) THEN
         IF(SFAC_S.GT.0.001) WRITE(6,95) NINT(ECHAR_S),
     >      TENFLX_S/6.25E11,SFAC_S
         IF(SFAC_S.LE.0.001) WRITE(6,96) NINT(ECHAR_S),
     >      TENFLX_S/6.25E11,SFAC_S
      ENDIF
 95   FORMAT(5X,'Gaussian AURORA on: Characteristic E='
     > ,I8,' eV; Energy Flux=',F9.2,' ergs/cm2/s; SFAC=',F8.2) 
 96   FORMAT(5X,'Maxwellian AURORA on: Characteristic Energy='
     > ,I8,' eV; Energy Flux=',F9.2,' ergs/cm2/s; SFAC=',F8.2)
      RETURN
      END    
